Gravitational Radiation from the Coalescence of Binary 
Neutron Stars: Effects Due to the Equation of State, Spin, and 

Mass Ratio 

Xing Zhuge, Joan M. Centrella, and Stephen L. W. McMillan 

Department of Physics and Atmospheric Science, Drexel University, Philadelphia, PA 19104 

On 

y—i ■ 

Abstract 

o 

t— i ; 

We calculate the gravitational radiation produced by the coalescence of 

> 

ON ' inspiraling binary neutron stars in the Newtonian regime using 3-dimensional 

m ' 
o 
o 



cr 



numerical simulations. The stars are modeled as polytropes and start out 



' in the point-mass regime at wide separation. The hydrodynamic integration 

On" 

*q | is performed using smooth particle hydrodynamics (SPH) with Newtonian 



gravity, and the gravitational radiation is calculated using the quadrupole 
approximation. We have run a number of simulations varying the neutron star 



^ ' radii, equations of state, spins, and mass ratio. The resulting gravitational 

waveforms and spectra are rich in information about the hydrodynamics of 
coalescence, and show characteristic dependence on GM/Rc 2 , the equation of 
state, and the mass ratio. 
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I. INTRODUCTION 



Coalescing binary neutron stars are among the most promising sources of gravitational 
waves for detection by interferometers such as the Laser Interferometric Gravitational-wave 
Observatory (LIGO) 0, VIRGO @, and GEO 0. Recent studies [§ suggest that binary 
inspiral due to gravitational radiation reaction, and the eventual coalescence of the com- 
ponent stars, may be detectable by these instruments at a rate of several per year. The 
inspiral phase comprises the last several thousand binary orbits and covers the frequency 
range / ~ 10-1000Hz, where the broad-band interferometers are most sensitive. During 
this stage, the separation of the stars is much larger than their radii and the gravitational 
radiation can be calculated using post-Newtonian expansions in the point-mass limit [1. 
Analysis of the inspiral wave form is expected to reveal the masses and spins of the neutron 
stars, as well as the orbital parameters of the binary systems f|[7|||f| . 

When the binary separation is comparable to the neutron star radius, hydrodynamic 
effects become dominant and coalescence takes place within a few orbits. The coalescence 
regime probably lies at or beyond the upper end of the frequency range accessible to broad- 
band detectors, but it may be observed using specially designed narrow band interferometers 
pTjj] or resonant detectors [PJ| . Such observations may yield valuable information about 
neutron star radii, and thereby the nuclear equation of state [|,|T^,|T^]. 

Three-dimensional numerical simulations are needed to study the detailed hydrodynam- 



ical evolution of the coalescence. Shibata, Nakamura, & Oohara jl^JT5|1 have studied the 
behavior of binaries with both synchronously rotating and non-rotating stars, using an Eu- 



lerian code with gravitational radiation reaction included. Ruffert, et al. |16[ have also used 
Eulerian methods with radiation reaction included to study coalescence of neutron stars with 
a physical equation of state and various spins. Rasio & Shapiro |T7|JT8|] have simulated the 
coalescence of synchronously rotating neutron-star binaries using the Lagrangian smooth 



particle hydrodynamics (SPH) method with purely Newtonian gravity. Davies et al. [19| 
have carried out SPH simulations of the inspiral and coalescence of nonsynchronously rotat- 
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ing neutron stars, focusing on the thermodynamics and nuclear physics of the coalescence. 
All of these studies used the quadrupole approximation to calculate the gravitational radi- 



ation emitted. Finally, Wilson, Mathews, and Marronetti [2(J have developed an Eulerian 
code that incorporates general relativistic effects in the limit in which the metric remains 
conformally flat and gravitational radiation is neglected. A multipole expansion is used to 
calculate the gravitational radiation. 

We have carried out 3-D simulations of binary neutron star coalescence in the Newtonian 
regime using SPH, with particular application to the resulting gravitational wave energy 
spectrum dE/df . The neutron stars are initially modeled as spherical polytropes on circular 
orbits, with separations sufficiently large that tidal effects are negligible. The stars thus start 
out effectively in the point-mass regime. The gravitational field is purely Newtonian, with 
the gravitational radiation calculated using the quadrupole approximation. To cause the 
stars to spiral in, we mimic the effects of gravitational radiation reaction by introducing a 
frictional term into the equations of motion to remove orbital energy and angular momentum 
at the rate given by the equivalent point-mass inspiral. As the neutron stars get closer 
together the tidal distortions grow and eventually dominate, and coalescence quickly follows. 
The resulting gravitational wave forms match smoothly onto the point-mass inspiral wave 
forms, facilitating analysis in the frequency domain. In Paper I [21] we considered equal mass 



neutron stars with M = 1.4 M Q and varied the neutron star radius and equation of state. 
We demonstrated that the resulting gravitational wave signatures are rich in information 
about the hydrodynamics of coalescence and are sensitive to both GM/Rc 2 and the equation 
of state. In this paper, we extend our study to include the effects of unequal masses as well 
as spin. 

It is important to understand the context of these models. Our work is carried out in 
the Newtonian regime and therefore is a first step toward understanding the gravitational 
radiation signatures of binary coalescence. With these simplified models we are able to study 
binaries that start out with fairly wide separations and make ~ 3 orbits before contact, and 
to concentrate on the hydrodynamical properties of the merger. Of course, the Newtonian 



approximation does break down for systems involving neutron stars, since GM/Rc 2 ~ 0.2 
for a typical neutron star of mass M = 1.4M and radius R = 10km. General relativistic 
effects can therefore be expected to play an important role in the final stages of inspiral and 



coalescence [|22|| , and Newtonian results must be viewed with appropriate caution. We believe 
that Newtonian models provide an interesting first look at the properties of coalescence 
waveforms and spectra. In addition, they can be used for comparison with general relativistic 
calculations to help determine where relativistic effects become important and how they show 
up in the resulting gravitaional waveforms and spectra. Finally, the valuable experience 
gained in carrying out these Newtonian calculations is important for the development of 
fully general relativistic models. 

This paper is organized as follows. In Sec. H we present a brief description of the 
techniques used in our simulations. The use of frictional terms in the equations of motion 



to mimic the effects of gravitational radiation reaction is discussed in Sec. |TT|. Section [TV 



revisits the standard model (with identical neutron stars having masses M = 1.4M and radii 
R = 10km), extending and expanding the analysis begun in Paper I. The effects of changing 
the neutron star radius, equation of state, and spin are examined in Sec. |V|. Binaries with 
unequal mass components are considered in Sec. [VT|. Finally, Sec. |V11| contains a summary 
and discussion of our results. 

II. SIMULATION TECHNIQUES 

The methods we used to produce our models have been presented in some detail in Paper 



I |2T| and reference p3[ . We therefore give a only brief description of these methods in this 
section, and refer the reader to the literature for further information. 

Lagrangian techniques such as SPH |24| are especially attractive for modeling neutron 
star coalescence since the computational resources can be concentrated where the mass 
is located instead of being spread over a grid that is mostly empty. We have used the 



implementation of SPH by Hernquist & Katz [^3J known as TREESPH. In this code, the 
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gravitational field is purely Newtonian and a hierarchical tree method [26| optimized for 



vector computers is used to calculate the gravitational forces. This leads to a significant 
gain in efficiency and allows the use of larger numbers of particles than would be possible 
with methods that simply sum over all possible pairs of particles. 

We calculate the gravitational radiation quantities in the quadrupole approximation, 



which is valid for nearly Newtonian sources [27]]. The reduced (i.e., traceless) quadrupole 
moment of the source is given by 

hi = J P ~ \5ijr 2 ) d 3 r, (1) 

where i, j = 1, 2, 3 are spatial indices and r = (x 2 + y 2 + z 2 ) 1//2 is the distance to the source. 
For an observer situated on the axis at = 0,0 = of a spherical coordinate system with 
its origin located at the center of mass of the source, the gravitational wave amplitudes for 
the two polarization states are given by 

G 1 

h + = -^~(hx - lyy), (2) 

G 2 •■ 

hx = ~~7~Ixy (3) 

Here, an overdot indicates a time derivative d/dt. The standard definition of gravitational- 
wave luminosity is 

where there is an implied sum on i and j, the superscript (3) indicates the third time 
derivative, and the double angle brackets indicate an average over several wave periods. 
Since such averaging is not well-defined during coalescence, we simply display the unaveraged 
quantity (G/5c 5 )l[f l[f in the plots below. The gravitational wave energy spectrum dE/df, 
which gives the energy emitted as gravitational radiation per unit frequency interval, is a 
key diagnostic tool for understanding the results of our simulations. It is given by Thorne 
|28| in the form 

^ = ^|(4vrr 2 )/ 2 (|/ i+ (/)| 2 + |^(/)| 2 ), (5) 



where h(f) is the Fourier transform of h(t), and the angle brackets denote an average over 
all source angles. See Paper I for details. 



We use the techniques of [£3|] to calculate the reduced quadrupole moment fjj and its 
derivatives. In particular, and are obtained using particle positions, velocities, and 
accelerations already present in the code to produce very smooth wave forms. This yields 
expressions similar to those of Finn and Evans p9[. However, iff requires the derivative 
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of the particle accelerations, which must be determined numerically and introduces noise 
into the gravitational wave luminosity L. We have applied smoothing to reduce this noise 
in producing all graphs of L in this paper; see [^] for further discussion. 

The neutron stars are initially modeled as widely separated polytropes with equation of 
state 

P = Kp ? = Kp 1+l / n , (6) 

where K is a constant that measures the specific entropy of the material and n is the 
polytropic index. The stars are placed on orbits with wide enough separation that tidal 
effects are negligible. An individual star may be allowed to be in uniform rotation about 
an axis through its center of mass. We take the direction of this spin angular velocity fl s 
(which is measured in an inertial frame) to be either parallel or anti-parallel to the direction 
of the orbital angular momentum. Because the time scale for tidal effects to develop is far 
greater than the dynamical time to for an individual star, where 

we start with stable, "cold" polytropes. The nonrotating stars (Q s = 0) were produced by 
the method discussed in |23|] . The rotating stars were produced using the method described 



in |30] and |5T . 
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III. MODELING INSPIRAL BY GRAVITATIONAL RADIATION REACTION 



Widely separated binary neutron stars (that is, with separation a 3> R) spiral together 
due to the effects of energy loss by gravitational radiation reaction. Once the two stars 
are close enough for tidal distortions to be significant, hydrodynamical effects dominate 
and rapid inspiral and coalescence ensue. In our calculations the neutron stars are placed 
on (nearly) circular orbits with wide enough separation that the stars are effectively in 
the point-mass limit. Since the gravitational field is purely Newtonian and does not take 
radiation reaction into account, we must explicitly include these losses to cause inspiral until 
hydrodynamical effects take over. 

We accomplish this by adding a frictional term to the particle acceleration equations 



to remove orbital given by the point-mass inspiral expression (see |19| for 

a similar approach). The gravitational wave luminosity for point-mass inspiral on circular 
orbits is El 



dE 
pm ~ ~dt 



32 G 4 ^M* {s) 



5 c 5 a 5 

pm 



where M. = Mi + M 2 is the total mass of the system, \i = M\M 2 j ' M. is the reduced mass, 
and the subscript "pm" refers to point-mass inspiral. We assume that this energy change 
is due to a frictional force / that is applied at the center of mass of each star, so that each 
point in the star feels the same frictional deceleration. For star 1, we obtain 



Mt \ - 1 dE 
~dt 



(9) 



pm 



where V\ is the center of mass velocity of star 1; an analogous expression in which the 
subscripts "1" and "2" are interchanged holds for star 2. Since f\ acts in the direction 
opposite to V\ this gives an acceleration 



(10) 



pm 

for star 1 and similarly for star 2. These frictional terms are added to the particle acceleration 
so that all particles in a given star experience the same frictional deceleration. The net 



effect is that the centers of mass of the stars follow trajectories that approximate point-mass 
inspiral. The frictional terms are applied until tidal effects dominate, leading to more rapid 



inspiral and coalescence by purely Newtonian hydro dynamical processes ||33|| . For each of 
the runs reported in this paper, we determine the optimal time to turn off the frictional 
terms experimentally; see Paper I for details. (Operationally, our assignment of a particle 
to a "star" is based simply on which body it happened to belong to initially. Since the 
frictional term is turned off before coalescence occurs, the question of what to do after the 
stars have merged does not arise.) 

The stars are initially placed on the x axis on a counter-clockwise circular orbit with 
separation a® in the center of mass frame of the system in the x — y plane. Thus, the center 
of mass of M\ is located at (x, y) position (a 1; 0) and that of M 2 is located at {—a 2 , 0), where 
a = «i + a 2 , d\ = aofi/Mi and similarly for a 2 . fl32| . The stars are then given the equivalent 
point-mass circular velocities V yA = (G/Ma ) l / 2 M 2 and V V)2 = -(G/MaoY^M^ 

To ensure that the stars start out on the correct point-mass inspiral trajectories, we 
also give them an initial inward radial velocity V x as follows. For point-mass inspiral the 
separation a(t) is given by 

a(t) =a (l- -J , (11) 
where ao is the separation at the initial time t = and 

T ° = 2o6 h ^ (12) 
is the inspiral time, i.e. the time needed to reach separation a = 0. We write 

da 



J*%^. (13) 

Requiring the center of mass of the binary system to have zero velocity then gives V Xj i = 
— (M 2 / M)V T and V Xy2 = {M\J M)V r - The use of the correct initial inspiral trajectory allows 
us to match our gravitational wave forms smoothly to the equivalent point-mass wave forms. 
This is important when analyzing the signals in the frequency domain. 
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IV. BINARY COALESCENCE: THE STANDARD MODEL 



We begin by examining binary coalescence for the case of equal mass neutron stars with 
masses M 1 = M 2 = M = 1.4M , radii R x = R 2 = R = 10km (so GM/Rc 2 = 0.21), and 
polytropic index n = 1 (T = 2). The stars start out with zero spin (O s = 0) on (nearly) 
circular orbits with initial separation a = 40 km in the point-mass limit. Time is measured 
in units of the dynamical time tn for a single star using equation (|7|); here, to = 7.3 x 10~ 5 s. 
We consider this to be our standard run and refer to it as Run 1; the parameters of this model 
are summarized in Table |. This run was first discussed in Paper I, which also included tests 
of our numerical method with varying particle number and artificial viscosity coefficients. 
In this section we revisit the standard model, extending and expanding the analysis begun 
in Paper I, defining our terminology, and reintroducing some key features of the problem. 



A. General Features of the Coalescence 

The evolution of this model with N = 4096 particles per star is shown in Figure [1]. In 
each frame all particles are projected onto the x — y plane. As the stars spiral together, 
their tidal bulges grow. By t = IOO^d, the center-of-mass separation of the two stars is 
~ 2.5R. At this point the stars undergo a dynamical instability driven by Newtonian tidal 
forces [J33|, causing the stars to fall together faster than they would on point-mass orbits. 



We therefore turn off the frictional term in the code at t — IOO^d and follow the rest of the 
evolution using purely Newtonian hydrodynamics and gravity. The stars rapidly merge and 
coalesce into a rotating barlike structure. Spiral arms form as mass is shed from the ends 
of the bar. Angular momentum is transported outward by gravitational torques and lost to 
the spiral arms. The arms expand and merge to form a disk around the central object. At 
the end of the run, the system is roughly axisymmetric. 

The tidal interactions between the stars increase as they spiral together. Even in the 
absence of fluid viscosity, Lai and Shapiro |34| have shown that a dynamical tidal lag angle 
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adyn develops due to the finite time needed for the structure of the stars to adjust to the 
rapidly changing tidal potential. This leads to the formation of tidal bulges that are not 
directly aligned, and the resulting gravitational torques cause each star to spin. As the 
stars spiral together and the separation decreases, ad yn becomes larger. In Figure ^ (c) we 
estimate the lag angle of the stars near contact to be ~ 15°. For non-spinning stars with the 
parameters of Run 1, Lai and Shapiro |34| find «d y n ~ 12° at contact, in good agreement 



with our results. 

Now consider the coalescence in a reference frame co-rotating with the binary. As seen 
from this co-rotating frame, the stars appear to be spinning in opposite directions. As the 
stars begin to merge, the fluid in star 1 is moving in the direction opposite to the fluid in star 
2 at the point of contact. This velocity difference over a short spatial scale gives rise to a shear 
layer between the stars, which is subject to the Kelvin-Helmholtz instability p5| , p6| . For 



real neutron stars, vortices can form within this layer on small scales; these eddies can grow 
and merge together, and turbulence can develop. The behavior of this turbulent region can 
be important in determining the mixing of the material from the individual neutron stars to 
produce a final remnant. For example, a turbulent viscosity could be generated and thereby 
the final configuration may not be irrotational. 

However, these processes are difficult to model accurately using 3-D numerical simula- 



tions due to their limited resolution |37],[^| . For example, vortices can form on spatial scales 
determined by the resolution of the model and spurious numerical shear viscosity, rather 
than on the smaller physical scales expected in real fluids. These shortcomings must be 
taken into account when interpreting the results of numerical models of binary coalescence. 

We have examined our simulations to see if these effects are occurring. We find that, as 
the stars merge, a couple of eddies form across the shear layer where the two stars meet. 
By t ~ 127tD (Figure [I] [g]), the stars have coalesced to the point that there is only a 
single density maximum at the center of the remnant. Each star in Run 1 has N = 4096, 
giving an effective resolution of ~ N 1 ^ 3 = 16 particles across the diameter of each star. 
The Eulerian models of Ruffert, Janka, and Schafer 0] have somewhat better resolution 
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with the diameter of a single star covering ~ 20 — 40 zones. Their calculations show the 
development of two eddies across the shear layer. Recent work by Rasio P5| , using SPH with 
a larger N and increased resolution of the shear layer due to the placement of many smaller 
mass particles in the outer regions of each star, shows the development of more eddies on 
smaller scales. 

For these reasons, we suspect that numerical effects may be influencing the behavior 
of the shear layer in our models. For example, they might cause the stars to mix more 
rapidly than they would in reality. Also, the bar phase of the evolution might be of shorter 



duration, and the final remnant might have different properties. However, Ruffert, et al. fl6 
point out that, in real neutron stars, the Kelvin- Helmholtz mechanism may develop into the 
macroscopic regime on the spatial scale of the coalescing stars as the eddies grow and merge. 
This could produce a final flow pattern similar to that seen in the simulations. They also 
suggest alternative explanations for the development of the final flow pattern. Clearly, more 
work is needed to resolve these issues. In particular, simulations with many more particles 
(and more grid zones in the Eulerian case) need to be done. 

The two neutron stars merge to form an object of total mass ~ 2.8 M , with ~ 94% 
(~ 2.6 M ) of the matter and ~ 74% of the angular momentum in a central core w < 2R, 
and the remainder in a disk (see Paper I). As the evolution proceeds, the core mass and 
angular momentum remain nearly constant, while the angular momentum in the disk is 
transported outward due to the effects of gravitational torques. By the end of the simulation 
at t = 200£d, the central core is axisymmetric with material on the edge of the remnant at 
V3 ~ 2R rotating at Q ~ 0.6^, near the critical angular velocity for breakup. We believe the 
reason for this axisymmetry can be understood as follows. Consider an equilibrium sequence 
of uniformly rotating axisymmetric polytropes parametrized by (3 = T TOt /\W\, where T rot is 
the rotational kinetic energy and W is the gravitational potential energy. As j3 is increased 
along such a sequence, a point is eventually reached at which mass is lost at the equator. 
Uniformly rotating polytropes with n > 0.808 (r < 2.24) reach this mass-shedding limit 



before the point at which ellipsoidal configurations can exist |39|;|4(|. Although the central 
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core in our simulation does show differential rotation, we believe that a similar mechanism 
is operating here, causing the core to be essentially axisymmetric at the end of the run. 

Stable, nonrotating neutron stars are believed to have a maximum mass in the range 
~ 1.4 M to ~ 2.2 M Q pl| , ^2| , depending on the equation of state. Rotation can increase 
this by up to ~ 17%, yielding a maximum mass < 2.6 M |43| , |44| (again depending on the 
equation of state), for rotation near breakup speed. Because the gravitational field in these 
models is purely Newtonian, the merged remnant in these simulations cannot collapse to a 
black hole. However, since the final core mass is greater than or comparable to the maximum 
allowed neutron star mass, general relativity may cause a black hole to form. 

The ability of rotation to prevent gravitational collapse to a black hole can be estimated 
by examining the dimensionless parameter A = cJ/GAi 2 , where Ai refers to the mass of 



the entire system. Piran and Stark |45| modeled the gravitational collapse of rigidly rotating 
polytropes with T = 2 and < A < 1.5 using a fully general relativistic 2-D axisymmetric 
code. They found that for A < A cr i t , the collapsing object formed a black hole, whereas for 
A > A cr it the collapse was halted by centrifugal forces leading to a bounce with no black 
hole formation. For the cases they considered, they found ^. cri t ~ 0.8 — 1.2. 

It is interesting to estimate the value of A for our simulations. Since this rotation 
parameter is a general relativistic concept and the values of J and M. that we determine 
from our simulations are purely Newtonian, this discussion must be treated with caution. 
With these caveats, we can examine Figure 0, which shows A versus cylindrical radius w 
for several times during the coalescence of Run 1. At t — 127tD, the innermost regions 
have A > 1. As the evolution proceeds, gravitational torques transport angular momentum 
outward and the central value of A drops. In all cases, A < A C iit at the core radius w ~ 2R. 
Thus, our Newtonian results indicate that it is likely a black hole will form. Of course, a 
firm determination of the final result of binary neutron star coalescence must await a fully 
general relativistic calculation. 
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B. Properties of the Emitted Gravitational Radiation 



Figure |3](a) shows the gravitational wave form rh + for an observer located on the axis at 
6 = = at distance r from the source. (For simplicity, we show only one polarization, rh + . 
For all runs presented in this paper, rh x is very similar in appearance to rh + , with a phase 
shift of 90°.) The solid line gives the code wave form and the dashed line the point-mass 
result. The code wave form matches the point-mass case for the first couple of orbits (note 
that the orbital period is twice the gravitational wave period, T or bit = 2Tgw)- As the tidal 
bulges grow and the stars spiral together faster than they would on point-mass trajectories, 
the gravitational waves increase in both amplitude and frequency (cf. p3|). Figure ||(b) 
gives the gravitational wave luminosity L/Lq, where Lq = c 5 /G. The code results (solid 
lines) initially track the point-mass case (dashed lines), then depart significantly from the 
point-mass predictions somewhat before the onset of dynamical instability. The amplitudes 
of the wave forms and luminosity both reach their maximum values during the early stages 
of the merger at t ~ 105 — IIO^d, when the numerical effects discussed above are least 
important. The amplitudes then decrease as the coalescence proceeds. By t ~ 180£d, the 
gravitational waves have shut off and the system is essentially axisymmetric. 

Table [TI] gives the maximum amplitudes of the gravitational wave forms and luminosities 
for the runs presented in this paper. For Run 1, we find that the maximum value of the wave 
form for a source located at distance r from the observer is (c 2 /GM)r |/i max | ~ 2.0(GM/Rc 2 ). 
For comparison, Rasio and Shapiro ( ||18|| ; this is listed as entry RSa in our Table [n|) 
found (c 2 /GM)r |/i max | ~ 2A(GM/Rc 2 ) for a synchronous binary with T = 2. Also, 
the maximum luminosity is (L max /L ) ~ 0.39(GM/Rc 2 ) 5 . Rasio and Shapiro found 
(L max /L ) ~ 0.55(GM/Rc 2 ) 5 . 

The gravitational wave energy spectrum dE/df has proved very useful for analyzing the 
models in the frequency domain. For point-mass inspiral dE/df ~ f^ 1 ^ 3 [p8fl , where the 
decrease in energy with frequency is due to the fact that the binary spends fewer cycles 
(and hence emits less energy) near a given frequency as it spirals in. Although our runs do 
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start out in the point-mass regime, the stars merge and coalesce within just a few orbits. To 
achieve a reasonably long region of point-mass inspiral in the frequency domain, we match 
the code wave forms for all runs in this paper onto point-mass wave forms extending back 
to much larger binary separations and thus lower frequencies. See Paper I for details. 

Figure |3](c) shows the energy spectrum dE/df for Run 1. The short dashed line is the 
spectrum for the wave forms truncated at time t = 120£d; this probes the initial stages of 
the merger, during which any numerical effects due to Kelvin-Helmholtz mixing should be 
less important. The solid line shows the spectrum for the full (i.e. untruncated) wave forms. 
In Paper I we defined a number of characteristic frequencies based on features observed 
in the energy spectrum, and identified these features with various dynamical phases of the 
coalescence. Starting in the point-mass regime, as / increases, dE/df first drops below the 
point-mass inspiral value and reaches a local minimum at / ~ 1500Hz. Recall that, for 
point masses at separation a, the gravitational wave frequency (which is twice the orbital 
frequency) is given by 

, 1 / GM \ 1/2 n .. 

/cw, pm = - (—) • (14) 

For point-mass inspiral, the frequency at separation a ~ 2.5R (where dynamical instability 
is predicted to set in) is fd yn ~ 1500Hz. We therefore identify this dip with the onset of 
dynamical instability. 

Beyond /d yn j the spectrum for the truncated wave forms (short dashed line) rises to a 
local maximum at / ~ 2000 Hz, then drops off sharply at higher frequencies. Using equation 
(fUj), we see that the point-mass gravitational wave frequency at "contact" (i.e. at separation 
a = 2R) is /contact ~ 2200Hz. We associate the peak with the formation of a rapidly rotating 
merged system. 

The solid line in Figure 0(c), which shows dE/df for the entire wave forms, contains 
features characterizing the later stages of the coalescence. We see that the peak shifts 
to higher frequencies, / pea k ~ 2500Hz, which we attribute to the formation of a transient, 
rotating barlike structure as the stars coalesce, between t = 120£d and t = 150tD- Continued 
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shrinking of the merged system as the coalescence proceeds causes the rotation speed of 
the bar to increase. We note that, if numerical effects due to Kelvin-Helmholtz mixing 
are shortening the duration of the bar phase of the evolution, our simulations could be 
substantially underestimating the strength of this feature. 

Beyond f pea k, the (untruncated) spectrum drops sharply, then rises to a secondary max- 
imum at / S ec ~ 3200Hz. It appears that this feature is due to transient oscillations induced 
in the merged remnant during coalescence, and therefore this region of the spectrum is un- 
doubtedly influenced by the numerical effects discussed above. We point out that, whether 
or not the final outcome of actual neutron star coalescence is a black hole, it seems reasonable 
that a spectrum of quasinormal mode oscillations will be produced as the remnant "rings 
down" . Higher resolution simulations that are able to track the detailed mixing of the stars 
during coalescence, and that fully incorporate general relativity, are needed to calculate this 
region of the spectrum accurately. 

V. BINARIES WITH EQUAL-MASS COMPONENTS 

Gravitational radiation from the coalescence of binary neutron stars is expected to con- 
tain important information about both the stars themselves and the interaction between 
them. We have parametrized the binary components as polytropes of masses Mi and M2, 
radii R\ and R2, spins Q 8j i and fl St 2, and equation of state T. In this section we present the 
results of runs with components having equal masses Mi = M 2 = M — 1.4M and radii 
By = i? 2 = R, and varying R, T, and Q s . The parameters of these runs are summarized in 
Table |. 

A. Varying the Neutron Star Radius R 

Run 2 is the same as Run 1 except that the neutron star radius R = 15km; with 
M = 1.4M , this gives GM/Rc 2 = 0.14. Aspects of this run have already been presented 
in Paper I; we include it here for completeness. In Run 2, the stars start out at initial 
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separation ao = 4i2 = 60km and the gravitational friction terms are turned off at t = 245tD? 
just before the stars come into contact. As in the standard run rapid coalescence takes place, 
with mass shed through spiral arms, producing a roughly axisymmetric final remnant. 

Figure f|(a) shows the gravitational wave form rh+ for an observer on the axis at 9 — <fi — 
at distance r from the source. The gravitational wave luminosity L is shown in Figure ^(b). 
The maximum amplitudes of the wave forms and luminosity both occur at t ~ 260£d- The 
maximum value of the wave form (see Table |T|) for a source located at distance r from the 
observer is (c 2 /GM)r\h 

max| °° 2. 1 {GM/Rc 2 ^ and the maximum luminosity is (Z/ max /Zvo) ~ 
0. 39 (GM/Rc 2 ) 5 ; both of these peak amplitudes are essentially the same as those obtained 
in the standard run. 

The gravitational wave energy spectrum dE/df for Run 2 is shown in Figure ^(c). For this 
run, dynamical instability is expected to occur at separation a ~ 2.5R [[33], corresponding to 
a point-mass inspiral frequency fd y n ~ 850Hz. As before, this behavior is evident in the data 
as the spectrum drops below the point-mass value, reaching a minimum at ~ 800 — 900Hz. 
The truncated spectrum (dashed line) then rises to a local maximum at / ~ 980Hz before 
dropping off at higher frequencies. For comparison, the point-mass frequency at separation 
a ~ 2R is /contact ~ 1200Hz. The solid line in Figure (|(c) shows dE/df for the full run; once 
again, the peak becomes more pronounced and moves to higher frequencies (/peak ~ 1350Hz) 
when the late-time wave form is included. Note that the spectrum does not rise above the 
point-mass result as in Run 1. Nevertheless, it drops sharply just beyond / pea k, rising again 
to a secondary maximum at / sec ~ 1800Hz. The lack of a strong peak may be due to the 
weaker tidal forces at the point of dynamical instability, which occurs at a larger physical 
separation than in Run 1, producing a less-pronounced and shorter-lived bar. As discussed 
above, numerical effects may well also be weakening this feature. 

It is quite instructive to compare the spectra of Runs 1 and 2. Since equation (|TJ]) gives 
the frequency of gravitational radiation for point-mass inspiral ~ a~ 3 / 2 ~ R~ 3 ^ 2 at contact, 
we expect that the characteristic frequencies in the Run 2 spectrum should scale relative to 
those in Run 1 roughly as /2//1 ~ 0.54. Our results for / pea k and / sec do show this behavior. 
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B. Varying the Equation of State T 



Run 3 is the same as Run 1, except that we use a stiffer polytropic equation of state 
r = 3 (n = 0.5). A low-resolution version of this run (with iV = 1024 particles per star) 
was originally presented in Paper I; here we repeat this model with N = 4096. 

The key difference between the coalescence of equal-mass neutron stars in Runs 1 and 3 is 
that the core of the final remnant in Run 3 is non-axisymmetric. This result was also found 
by Rasio and Shapiro |T8|] for the case of synchronous binaries with stiff equations of state. 
As remarked above, uniformly rotating polytropes with T > 2.24 can sustain ellipsoidal 
shapes, since the mass-shedding limit along a sequence of equilibrium models is reached 
after the point at which the ellipsoidal sequence bifurcates from the spheroidal sequence. 
Although the coalesced remnant in Run 3 does have some differential rotation, we expect 
that this mechanism is operating in this model. In particular, the movement of mass out 
of the central core region through spiral arms indicates that the system is rotating at the 
mass-shedding limit. This process is shown in Figure p] (e) - (i) for Run 1. In Run 3, the 
spiral arms are somewhat narrower, as expected for a stiffer equation of state, and the core 
of the merged remnant is slightly non-axisymmetric. 

The signature of this rotating, non-axisymmetric core is clearly seen in the gravitational 
wave form rh + , shown in Figure [5](a). After reaching a maximum value at t ~ 160tD, 
the amplitude initially drops rather sharply, and then decays more gradually on a much 
longer timescale. The luminosity L also reaches its maximum value at t ~ 160to, and 
decreases slowly at late times (see Figure ||[b]); we find (c 2 /GM)r \h max \ ~ 1.9(GM/Rc 2 ) 
and (L max /L ) ~ 0.29(GM/Rc 2 ) 5 (see Table 0). 



For this run, dynamical instability is expected to occur at separation a ~ 2.76R 33 



corresponding to a point-mass inspiral frequency fd yn ~ 1340Hz. Figure ||(c) shows that 
the spectrum dE/df for the early stages of the merger (dashed line) drops below the point- 
mass value and reaches a local minimum near /dyn- It then rises to a local maximum at 
/ ~ 1800Hz, reaching an amplitude above the point-mass value. The spectrum for the entire 
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run (solid line) has a very pronounced peak at / pea k ~ 2200Hz due to contributions from 
the rotating, non-axisymmetric core. 

As T increases the polytrope becomes less centrally condensed, with a larger fraction of 
its mass in its outer regions. Therefore, tidal effects between the two stars become important 
at larger separations for T = 3 than for T = 2 [l"B| . Dynamical instability is thus predicted 



to occur at a larger physical separation in Run 3 than in Run 1, and we expect that the 
spectral features will appear at correspondingly lower frequencies. Using point-mass inspiral 
we estimate the ratio to be /3//1 ~ 0.86, and our results do indeed show this behavior. 

Run 4 is the same as Run 1, except that T = 5/3 (n = 1.5) and we used N = 2048. 
We note in passing that the equation of state for nuclear matter is generally believed to 
be stiff, with r > 2, and that the case T = 5/3 applied to compact objects may be more 
appropriate for low-mass white dwarfs |47[]. Nevertheless, this run allows us to explore the 
effects of a softer equation of state on binary neutron star coalescence, and we include it 
here for comparison. 

Rasio and Shapiro |I7J] have shown that, for a synchronous binary composed of equal- 
mass stars with T = 5/3, dynamical instability occurs at separation a ~ 2AR. At this time, 
the stars in their model were already in contact. In our case, the separation of the centers 
of mass of the two stars does drop sharply near a = 2AR, and the stars are just on the 
verge of contact at the time. As in Run 1, the stars in Run 4 merge to produce a rotating, 
axisymmetric final remnant. During this process, spiral arms form and move mass out of the 
central region into a flattened halo surrounding the central core. The spiral arms in Run 4 
appear somewhat wider than in Run 1, as expected, since a softer equation of state is used. 

The gravitational wave form rh + and the luminosity L are shown in Figure |^(a) 
and (b), respectively. Both quantities reach their maximum values at t ~ 115tD ; with 
(c 2 /GM)r\h ma3[ \ ~ 2.3(GM/Rc 2 ) and (L max /L ) ~ 0.59(GM/i?c 2 ) 5 . 

The spectrum dE/df for Run 4 is shown in Figure f](c). Comparison with Figure 0(c) 
shows that the spectrum in Run 4 follows the point-mass inspiral result to somewhat higher 
frequencies, and therefore smaller separations, than in Run 1. This is due to the fact that 



polytropes with T = 5/3 are more centrally condensed than those with T = 2, and hence 
approximate point masses to smaller separations. Using the value a = 2AR as the separation 
at which dynamical instability takes place, we find /dyn ~ 1650Hz. The spectrum does dip 
below the point-mass result to reach a shallow minimum around /dyn, then rises above it to 
a broad peak at / ~ 2100 — 2500Hz. We estimate / pca k for the full spectrum (solid line) to 
be /peak ~ 2300Hz. This is followed by a steep drop and a rise to a secondary maximum at 
/see ~ 4000Hz. If we again use the point-mass scaling / ~ a~ 3//2 , along with the above value 
of /dyn, we find /4//1 ~ 1.05. Our numerical results give ratios in the range ~ 1.05 — 1.20. 



C. Varying the Neutron Star Spin 0, s 

Thus far, we have studied neutron binary systems in which the components have zero 
spin at large separations (a > 4i?), as seen from a non- rotating frame. Of course, neutron 
star binaries are not expected to be synchronously rotating, since the time scale for syn- 
chronization is in general much longer than the time scale for orbital decay and inspiral 



due to gravitational radiation reaction ||48|| . However, it is very likely on both theoretical 
and observational grounds that neutron stars are born with non-zero spin, as evidenced by 
the many short-period pulsars known in the Milky Way galaxy. We have therefore run two 
models to investigate the effects of spin on the coalescence and the resulting gravitational 
wave signals. 

The stars used in these models are taken to be uniformly rotating with spin angular 
velocity Q s . They were produced using the method described in |3(J (see also |j3~i~]), with 
N = 4055. (In this method, one does not have complete control over the number of particles 
accepted into the star; this accounts for the somewhat unusual value of N.) We have 
chosen |fi s | = 0.175t5\ which is about 30% of the maximum rotation rate that a uniformly 
rotating neutron star can have before it sheds mass at its equator, |fi s max | < OM^, 1 P3| , ^9| . 
For M = 1.4M and R = 10km, this corresponds to a spin period T s ~ 2.6ms. The stars in 
our runs have negligible rotational flattening, with the polar radius ~ 98% of the equatorial 
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radius. Their spins are allowed to be either positive (i.e. parallel to the orbital angular 
momentum) or negative. 

In Run 5, both stars are spinning in the positive sense, with Q St i = fi Sj 2 = 0.175t5 1 at 
initial separation ao = 4R. The wave form rh + reaches its maximum amplitude at t ~ 105iD 
as shown in Figure 0(a). It attains a higher maximum amplitude than in Run 1, then drops 
abruptly to a considerably lower amplitude; cf. Figure ||](a). The luminosity, which is 
presented in Figure 0(b), also reaches a larger maximum value than in Run 1 at t ~ 105£d; 
see Table 0. We find (c 2 /GM)r\h max \ ~ 2.3(GM/Rc 2 ) and (L max /L ) ~ 0.59(GM/Rc 2 ) 5 . 
The spectrum dE/df is given in Figure 0(c). Overall, the spectral features are similar in 
appearance to those shown in Figure 0(c) for Run 1. However, both / pea k and / sec occur at 



somewhat higher frequencies in Run 5; see Table fU . 

Run 6 is the same as Run 5, except that the stars have spins f2 Sj i = — f2 Sj 2 = O.lTSt^ 1 . 
The gravitational wave form (Figure ||[a]) reaches its maximum amplitude at t ~ IOO^d- 
This maximum value is not as large as that obtained in Run 5, and the subsequent drop 
in the amplitude of the wave form is more gradual. The luminosity, shown in Figure |^(b), 
also reaches a maximum at t ~ lOOtn- Here, the maximum luminosity is also somewhat less 
than in Run 5, although greater than in Run 1, and there is no secondary maximum. In 
this case we find (c 2 /GM)r \h max \ ~ 2.0(GM/Rc 2 ), (L max /L ) ~ 0.39(GM/ Rc 2 ) 5 . Finally, 
Figure H(c) shows the spectrum dE/df. Comparison with Figure 0(c) shows that the spectral 
features for Run 6 are again similar in appearance to those of Run 1 and occur at about the 



same frequencies. These results are summarized in Table D ( 



VI. BINARIES WITH UNEQUAL-MASS COMPONENTS 

In this section we present the results of simulations with binary components of unequal 
masses. The primary is taken to have mass Mi = 1.4M and radius Ri = 10km. We 
then let the secondary have mass M2 = qM%, where q = M2/M1 < 1 is the mass ratio. 
To calculate the radius of the secondary R 2 , we follow Rasio and Shapiro [n| and assume 



20 



that the system has constant specific entropy. Thus, the initial state of both components is 
constructed using equation (0) with the same polytropic constant K x = K 2 and the same 
value of T. This gives the mass-radius relation 



- 



(15) 



Both components are taken to have zero spin, Q sA = f2 s 2 = 0. We use N = 4096, with all 
particles in a given star having the same mass. The dynamical time to used to describe the 
evolution of these runs is calculated using the parameters of the primary. Table | summarizes 
the parameters of these models. 

For this paper, we have carried out runs with two different mass ratios, q = 0.85 and q = 
0.5, using both T = 2 and T = 3. The value q = 0.85 is believed to be the most probable mass 



ratio for the binary pulsar PSR 2303+46 |50 |. Although this is the smallest observed value 



of q for a binary pulsar with a neutron star companion |1| , we consider it important at this 
early stage in our understanding of the gravitational wave signals from binary coalescence to 
explore more extreme mass ratios. For this reason we have also considered the case q = 0.5. 

Run 7 is a slightly asymmetric binary with q = 0.85 and T = 2. By equation (|T5D, both 
components have equal radii, R 2 = R\. In Figure |S] all particles are projected onto the 
x — y plane to show the evolution of this model. Tidal bulges develop as the stars spiral 
together, and mass transfer begins by t ~ 120tD- As the merger proceeds, a single spiral 
arm or elongated "tail" is formed from the secondary, seen in Figure [| (d) - (g). By the end 
of the simulation at t = 250tD, the secondary has completely merged with the primary. The 
rotating, axisymmetric remnant has a central core of radius ~ 1.5i?i and is surrounded by 
a flattened halo consisting of material from the secondary. 

The gravitational wave form rh + for Run 7 is shown in Figure |T0"|(a), and the luminosity 
L is shown in Figure |TU](t») . Both quantities reach their maximum amplitudes at t ~ 125^ 
during the early stages of the merger. The maximum value of the wave form (see Table |T1[) 
for a source located at distance r from the observer is (c 2 /GM)r |/i max | ~ 1.6(GM 1 /i? 1 c 2 ) 
and the maximum luminosity is (L max /L ) ~ 0.18(GMi/i?ic 2 ) 5 . Figure |l0|(c) shows the 
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gravitational wave energy spectrum dE/df for this run. 

In Run 8 q = 0.85 and T = 3 so that, according to equation R% = 9.7km. Figure |TT 
shows the evolution of this run, which proceeds in a similar way to that in Run 7. In this 
case, however, the stiffer equation of state leads to a narrower extended tail. The central 
region also retains an elongated shape for a longer time. By the end of the simulation at 
t = 200tD, the central core is essentially axisymmetric with radius ~ 1.5i?i and a flattened 
halo composed of mass from the secondary. The remnant is also slightly displaced from the 
sytem center of mass. Rasio and Shapiro also notice this behavior in their synchronous 
merger with q = 0.85 and T = 3, attributing it to the asymmetric, single-arm mass outflow. 
Interestingly, Run 7 with T = 2 also has asymmetric outflow, yet the final remnant suffers 
a much smaller displacement from the system center of mass. 

The gravitational wave form rh + for Run 8 is shown in Figure |T2"|(a), and the luminosity 
L is shown in Figure |12|(b). Both quantities reach their maximum amplitudes at t ~ IIO^d, 
during the early stages of the merger. More gravitational radiation is produced after the 
maximum is reached than in Run 7, due to the more strongly non-axisymmetric central 
region. The peak value of the wave form (see Table [TJ) for a source located at distance r from 
the observer is (c 2 /GM)r\h max \ ~ 1.6(GM 1 /fi 1 c 2 ); the maximum luminosity is (L max /L ) ~ 
0.19(GMi/i?ic 2 ) 5 . Figure |l2](c) shows the gravitational wave energy spectrum dE/df for 
this run. Comparing the spectrum for the full wave forms with that given in Figure |T0l(c) 
for Run 7, we see that Run 8 with T = 3 shows a more pronounced peak at a slightly lower 
frequency / pea k- This behavior is similar to that seen in comparing Runs 1 and 3 and is 



attributed to the stronger bar like central region; cf. § |VB . 

Run 9 has q = 0.5 and T = 2 which implies that the components have equal radii, 
i?2 = Ri- The evolution of this model is shown in Figure [13], where all particles are projected 
onto the x — y plane. As the stars spiral together, the secondary develops a tidal bulge and 
starts to lose mass to the primary. Although this mass transfer completely disrupts the 
secondary, the primary is not strongly affected. At the end of the run t = 300£d> most of 
the mass of Mi is spread in a flattened halo of radius ~ 5Ri around Mi, with a single spiral 
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arm extending out to very large radii, up to ~ 40i?i. The halo and core contain ~ 96% of 
the total system mass, and ~ 70% of the total angular momentum (relative to the center of 
mass of the system). 

Figure 0(a) shows the gravitational wave form rh + for this run, and Figure |14|(b) 
shows the luminosity L. Both the wave form and the luminosity reach their peak val- 
ues at t ~ 195£d, during the early stages of mass transfer. The gravitational waves 
shut off less than 1.5 orbits later, indicating that the secondary is quickly disrupted. 
The maximum value of the wave form (see Table ||) for a source located at distance r 
from the observer is (c 2 /GM)r |/i max | ~ 0.67(GM\/ Ric 2 ) and the maximum luminosity is 
(L max /L ) ~ 0.011(GMi/ Ric 2 )^ . The gravitational wave energy spectrum dE/df for this 
run is shown in Figure [14|(c). Since the evolution of the system proceeds rapidly by mass 
transfer, there is little difference between these two curves, with / pea k ~ 900Hz; there is no 
discernible high-frequency secondary feature present at late times. 

Run 10 uses q = 0.5 and r = 3, so R2 = 8.7km. The evolution of this model is shown 
in Figure ITS]. As before, the tidal bulge on the secondary grows and the resulting mass 



transfer starts to disrupt it. The primary feels little effect as matter from the secondary 
spreads around it. However, this period of mass transfer ends < 1.5 orbits after it began 
when the secondary, now much reduced in size, moves out to a wider orbit. At the end 
of the simulation (t = 300£d), the primary contains ~ 86% of the total system mass and 
has a radius ~ 1.5R\, where R\ is the initial primary radius. The secondary has a mass 
~ 0.29M Q , about 42% of its original mass, and is in orbit about the primary at center of 
mass separation roughly ~ 4.5i?i. Although we ended the simulation at this point, in reality 
the inspiral will begin again due to gravitational radiation reaction. 

The gravitational wave form rh + is shown for this run in Figure |16](a), and the luminosity 
in Figure 0(b). Both the wave form and luminosity attain their maximum amplitudes in 
the early stages of mass transfer at t ~ 220£d- At late times, the wave form shows the 
signature of the final binary orbit. The maximum value of the wave form (see Table [TT]) for 
a source located at distance r from the observer is (c 2 /GM)r\h max \ ~ 0.76(GMi/i?ic 2 ) and 
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the peak luminosity is (L max /L ) ~ 0.02(GMi/ J Ric 2 ) 5 . 

The gravitational wave spectrum dE/df for Run 10 is given in Figure 0(c). We find 
/peak ~ 900Hz. The spectrum for the full run shows a signal in the region around ~ 
600Hz. Since the point-mass gravitational wave frequency for the end state of the system 
at separation a ~ 4.5-Ri is ~ 610Hz, we believe that this feature is due to the final binary 
orbit. In addition, a high-frequency feature / sec ~ 2300Hz appears at late times. Since the 
amplitude of this secondary feature is so much smaller than the main signal around /peak, 
we suspect that it may be a numerical artifact. 

VII. SUMMARY AND DISCUSSION 

We have used SPH to perform 3-D simulations of the coalescence of binary neutron stars 
with the goal of determining the gravitational radiation signals produced and understanding 
the information that can be extracted from the wave forms, luminosities, and spectra. The 
stars are initially modeled as spherical polytropes with masses M\ and M2, radii R\ and R2, 
spins fi Sj i and f2 Sj 2, and equation of state T, as summarized in Table |. At the start of each 
run, the stars are placed on (nearly) circular orbits with wide separations so that the binaries 
are effectively in the point-mass limit. The gravitational field is purely Newtonian, and 
the gravitational radiation quantities are calculated using the quadrupole approximation. 
Frictional terms are added to the equations of motion to mimic the effects of gravitational 
radiation reaction and to cause the stars to spiral together. As the stars get closer, tidal 
effects grow and eventually dominate. At this point, the frictional terms are turned off and 
the coalescence proceeds by purely Newtonian hydrodynamics. 

Our first set of runs features binaries having components with equal masses Mi = M2 
\ I = 1.4M and radii R\ = R2 = R. We varied the radius R, equation of state T, and 
spin Q s . In all of these runs, coalescence occurs rapidly once the dynamical stability limit is 
reached. The merging stars form a rotating barlike structure, and spiral arms are produced 
as mass is lost from the ends of the bar. Gravitational torques cause angular momentum 
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to be transported outward and lost to the spiral arms. The arms expand supersonically 
and merge to form a disk around the central object. For the cases that we studied, the 
rotating core of the final merged remnant is axisymmetric for T = 5/3 and T = 2, and 
non-axisymmetric for T = 3. 

The gravitational wave signals for these runs start out following the point-mass results. 
As the tidal effects grow stronger and the stars begin to spiral in faster than they would on 
point-mass trajectories, the amplitudes and frequencies of the wave forms and the amplitude 
of the luminosities increase faster than the point-mass results. These amplitudes reach their 
maximum values in the early stages of the coalescence, soon after the stars come into contact. 
Table [II] shows the scaling relationships for these maximum values for our runs as well as 
four runs for synchronous binaries by Rasio and Shapiro [JlS | that are labeled RSa, b, c, and 



d. For non-spinning stars with T = 2, we find that increasing the value of T decreases the 
strength of these maxima; cf. RSa and RSb. Allowing the stars to have identical spins in 
the same direction as the orbital angular momentum also increases the maximum values. In 
particular, our Run 5 gives very similar results to RSa. Interestingly, our Run 6, in which the 
stars have spins that are equal in magnitude but opposite in direction, produces maximum 
values that are essentially the same as Run 1 with non-spinning stars. 

Examination of the wave forms and luminosities after the maximum values are attained 
also reveals certain trends. First consider the effects of changing T, starting with Run 1 
which has T = 2. When we decrease this parameter to T = 5/3 in Run 4, the waves shut 
off more abruptly than in Run 1, as can be seen by examining Figures |3](a) and |6|(a). We 
believe that this is due to the tendency for more compressible fluids (i.e., those with smaller 
values of T) to reach the mass-shedding limit at smaller values of ft = T Iot /\W\. However, 
when this parameter is increased to T = 3 in Run 3 the wave form amplitude falls off more 
gradually at late times due to the rotating, slightly non-axisymmetric core; cf. Figures 0(a) 
and |5](a). We attribute this to the ability of polytropes with T > 2.4 to sustain ellipsoidal 
shapes since the mass-shedding limit along a sequence of equilibrium models is reached after 
the point at which the ellipsoidal sequence bifurcates from the spheroidal sequence. When 
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r = 2 and the stars have parallel spins in the direction of the orbital angular momentum 
as in Run 5, the amplitude of the wave form drops quite abruptly after the maximum is 
reached; this is shown in Figure 0(a). Similar behavior was seen by Rasio and Shapiro fI7 



and Shibata, et al. |14]]. In Run 6 the spins are anti-parallel, and Figure |8](a) shows that 
the fall-off in the wave form is less abrupt than in the case of parallel spins, but more rapid 
than in the case of zero spin. In addition, the luminosity drops sharply after reaching its 
maximum value. All of these runs show three peaks with successively decreasing amplitudes 
in the luminosity except for Run 6, which has only one peak. The synchronous runs of 
Rasio and Shapiro [p!7||l8| show a single peak in the luminosity for T = 5/3,2, and 3, with 
a secondary peak appearing for T = 10. 

The gravitational wave energy spectrum dE/df is a useful tool for analyzing the models. 
In the point-mass regime, the spectrum falls off according to dE/df ~ f^ 1 ^ 3 p8 |. When 



tidal effects become dominant and the dynamical stability limit is reached, the spectrum 
drops below the point mass curve and reaches a local minimum around the corresponding 
frequency. For the early stages of the merger, the spectrum then rises to a broad local 
maximum, then falls off rather quickly at higher frequencies. At later times, the peak 
sharpens and moves to higher frequencies, due to a transient, rotating bar-like structure 
that forms during the coalescence. The spectrum next drops off very sharply, and then rises 
to a secondary maximum. While we can explain this secondary peak in terms of oscillations 



in the final remnant, it is unclear how reliable the simulation is at late times. Table |(I 
shows that all runs with R = 10km and T = 2 give very similar values for / pca k. When 
R = 15km, these features occur at lower frequencies which scale as expected, roughly as 
~ Ch an gi n g the equation of state to T = 3 and using R = 10km produces a somewhat 

smaller decrease in these frequencies, which we attribute to the occurrence of the dynamical 
instability at a slightly larger orbital separation. 

We performed two simulations with non-equal mass stars having mass ratio q = 0.85. 
In Run 7 with T = 2, the two stars merge to form a remnant that has a central core of 
radius ~ 1.5i?i surrounded by a flattened halo formed out of matter from the secondary. 
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The behavior of Run 8 with T = 3 is similar, except that the central rotating object remains 
elongated for a longer period of time. This is reflected in the gravitational wave spectrum 
dE/df shown as the solid line in Figure |P2|(c), which shows a pronounced peak at / pea k ~ 
2200Hz. In contrast, the spectrum for Run 7, which is shown as the solid line in Figure |10|(c) , 
has a much weaker feature at / pea k ~ 2300Hz. 

We have also carried out two runs with non-equal mass components having mass ratio 
q = 0.5. In Run 9 with T = 2, the merger proceeds by mass transfer that completely disrupts 
the secondary. The gravitational waves shut off less than 1.5 orbits after the maximum 
amplitude is reached. When T = 3, as in Run 10, mass transfer also initially disrupts the 
secondary. However, less than 1.5 orbits later this process stops as the secondary moves out 
to a wider orbit. At the end of the run the secondary has ~ 42% of its original mass and 
orbits the primary at center of mass separation ~ A.5R\, where Ri is the initial radius of the 
primary. Further evolution of this system will proceed on a secular timescale by gravitational 
radiation reaction. Rasio and Shapiro [|1^] also report the formation of a detached binary 
for the case of an initially synchronous system with T = 3 and q = 0.5. Interestingly, their 
maximum amplitudes for the wave form and luminosity are very similar to ours (although 
their masses and final orbital separation are not); cf. Runs 10 and RSd in Table |TI1 

It is important to understand the influence of numerical effects on the results of our 
simulations. In Paper I we investigated the use of various artificial viscosity coefficients. 
We also showed that our results do not depend strongly on N, the number of particles per 



star, for iV = 1024, 2048 and 4096. However, as discussed in § |V| above, we suspect that 
numerical effects due to Kelvin-Helmholtz instabilities may be influencing the behavior of the 
shear layer that forms where the two merging stars meet. In particular, the bar phase of the 
evolution may be artificially shortened and the final remnant may have different properties. 
These numerical effects should become severe after the maximum gravitational radiation 
amplitudes are reached. Thus, we believe that the scaling of the maximum amplitudes, the 
shape of the spectrum dE/df for the early stages of the merger, and the approximate location 
of the frequency / pea k are reliable. However, the shape of the gravitational radiation signals 
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after the maximum values are reached, the strength of the spectral feature at / pea k, and the 
secondary feature at / sec may be affected somewhat by numerical processes, and we advise 
caution when interpreting these results. We remark that the simulations of non- synchronous 
binaries carried out by other groups are also expected to suffer from these problems; see § [TV| . 
More work is needed to clarify these issues, including simulations with higher resolution. 

The gravitational waveforms and spectra resulting from these Newtonian simulations 
contain much information about the hydrodynamics of the merger. Of course, general rela- 
tivity is likely to bring in other physical effects that need to be studied in order to understand 



the data expected from the detectors. For example, Lai and Wiseman |22[ have recently 
shown that the inclusion of certain general relativistic effects along with Newtonian tidal 
processes causes the neutron stars to begin their final plunge towards merger at larger orbital 
separations, and hence at lower frequencies, than in the purely Newtonian case. Such infor- 
mation is potentially very important for the detection of the gravitational wave signals from 
binary neutron star coalescence by LIGO and other detectors. We intend to incorporate full 
general relativity in our future work. 
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TABLES 



111 1 1 

Model 


Q 


Ri 


R2 


a 




rji 

T s ,i 


T s , 2 


n 


i 


iV 






(km) 


(km) 


(km) 


(ms) 


(ms) 


(ms) 








Run 1 


1 


10 


10 


40 


0.073 








1 


2 


4096 


Run 2 


1 


15 


15 


60 


0.13 








1 


2 


1024 


Run 3 


1 


10 


10 


45 


0.073 








1/2 


3 


4096 


Run 4 


1 


10 


10 


40 


0.073 








3/2 


r In 

5/3 


2048 


Run 5 
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10 


10 


40 


0.073 


2.6 


2.6 


1 


2 


4055 


Run 6 


1 


10 


10 


40 


0.073 


2.6 


-2.6 


1 


2 


4055 


Run 7 


0.85 


10 


10 


40 


0.073 








1 


2 


4096 


Run 8 


0.85 


10 


9.7 


40 


0.073 








1/2 


3 


4096 


Run 9 


0.5 


10 


10 


40 


0.073 








1 


2 


4096 


Run 10 


0.5 


10 


8.7 


40 


0.073 








1/2 


3 


4096 



TABLE I. Parameters of the models are given. For all runs, M\ = 1.4Mq. The mass ratio 
is q = M2/M1, the stellar radii are R± and R2, and the initial separation is ao- The dynamical 
time tjy is calculated using the parameters of star 1. T s gives the spin period of each star, with 
a positive (or negative) value denoting a spin in the same (or opposite) direction as the orbital 
angular momentum. The polytropic index n and the L refer to the equation of state. Each star 
contains iV SPH particles. 
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Model (cVGMOrl/wl {c 2 / GM^h^^a' 1 10 4 (L max /L ) (L max /L )a" 5 



Run 1 


0.43 


2.0 


1.6 


0.39 


Run 2 


0.29 


2.1 


0.21 


0.39 


Run 3 


0.40 


1.9 


1.2 


0.29 


Run 4 


0.48 


2.3 


2.4 


0.59 


Run 5 


0.48 


2.3 


2.2 


0.59 


Run 6 


0.42 


2.0 


1.6 


0.39 


Run 7 


0.33 


1.6 


0.75 


0.18 


Run 8 


0.33 


1.6 


0.76 


0.19 


Run 9 


0.14 


0.67 


0.044 


0.011 


Run 10 


0.16 


0.76 


0.083 


0.020 


RSa 




2.4 




0.55 


RSb 




2.2 




0.37 


RSc 




1.6 




0.14 


RSd 




0.8 




0.018 



TABLE II. The maximum amplitudes of the gravitational wave form and luminosity are given 
for each model. The value (c 2 /GMi)r|/i max | ~ 0.4 corresponds to an amplitude h ~ 1.4 x 10~ 21 
for a source at distance r = 20Mpc (the approximate distance to the Virgo cluster of galaxies). 
Scaled versions of these quantities, in terms of the parameter a = GM\j R\c 2 , are also presented. 
The entries RSa,b, c, and d refer to models of synchronous binaries run by Rasio and Shapiro and 
are taken from Table I in reference [|l^]. Model RSa has V = 2 and 9=1; compare this with our 
Runs 1, 2, 5 and 6. Model RSb has T = 3 and q = 1; this should be compared with our Run 3. 
Model RSc has T = 3 and q = 0.85 and Model RSd has T = 3 and q = 0.5; compare these with 
our Runs 8 and 10, respectively. 
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Model 


/peak (Hz) 


/sec (Hz) 


Run 1 


2500 


3200 


Run 2 


1350 


1750 


Run 3 


2200 


2600 


Run 4 


2700 


4000 


Run 5 


2700 


3500 


Run 6 


2500 


3200 


Run 7 


2300 


3000 


Run 8 


2200 


2600 


Run 9 


940 




Run 10 


1000 


2300 


TABLE III. Various frequencies of the models relating to the spectra dE/df are given. The 
frequencies / pea k and / sec refer to the spectra calculated from the full wave forms and and give the 
bar rotation speed and the oscillation frequency for the remnant. The truncated spectra peak at 
somewhat lower frequencies. (The value of / pea k given here for Run 2 is somewhat smaller than 
that used in Paper I; the difference is due to some arbitrariness in estimating the location of the 



"cliff" in Fig. |(c).) 
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FIGURES 

FIG. 1. Particle positions are shown projected onto the x — y plane for Run 1. Here, 
M = 1.4M , R = 10 km, T = 2, and i D = 0.073 ms. The initial separation a = 4R. The 
stars are orbiting in the counter-clockwise direction. The vertical axis in each frame is y/R and 
the horizontal axis is x/R. 

FIG. 2. The dimensionless parameter A = cJ/GM 2 , where hA refers to the mass of the entire 
system, is shown as a function of cylindrical radius wjR for several times during the coalescence 
of Run 1. Here, R = 10 km as in Figure |l]. 

FIG. 3. (a) The gravitational wave form r/i + for Run 1 is shown for an observer located on the 
axis at 9 = <p = at distance r from the source. The solid line is the code wave form and the dashed 
line is the point-mass result, (b) The gravitational wave luminosity L/Lq, where Lq = c 5 /G. The 
solid line is the code result, and the dashed line gives the point mass profile, (c) The gravitational 
wave energy spectrum dE/df. The code wave forms have been matched onto point-mass results to 
produce a long region of point-mass inspiral at low frequencies. The solid line shows the spectrum 
for the entire run, and the short dashed line shows the spectrum for the wave forms truncated at 
t = 120£d- The long dashed is the point mass spectrum dE/df ~ f 1 ^- 

FIG. 4. (a) The gravitational wave form r/i + for Run 2 is shown for an observer located on 
the axis at 9 = (ft = at distance r from the source, (b) The gravitational wave luminosity L/Lq, 
where Lq = c 5 /G. (c) The gravitational wave energy spectrum dE/df. The solid line shows the 
spectrum for the entire run, and the dashed line shows the spectrum for the wave forms truncated 
at t = 270t D - 

FIG. 5. (a) The gravitational wave form rh+ for Run 3 is shown for an observer located on 
the axis at 9 = (ft = at distance r from the source, (b) The gravitational wave luminosity L/Lq, 
where Lq = c 5 /G. (c) The gravitational wave energy spectrum dE/df. The solid line shows the 
spectrum for the entire run, and the dashed line shows the spectrum for the wave forms truncated 
at t = 175tD- 
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FIG. 6. (a) The gravitational wave form rh + for Run 4 is shown for an observer located on 
the axis at 6 = (f> = at distance r from the source, (b) The gravitational wave luminosity L/Lq, 
where Lq = c 5 /G. (c) The gravitational wave energy spectrum dE/df . The solid line shows the 
spectrum for the entire run, and the dashed line shows the spectrum for the wave forms truncated 
at t = 125t D - 

FIG. 7. (a) The gravitational wave form rh + for Run 5 is shown for an observer located on 
the axis at 6 = 4> = at distance r from the source, (b) The gravitational wave luminosity L/Lq, 
where Lq = c 5 /G. (c) The gravitational wave energy spectrum dE/df. The solid line shows the 
spectrum for the entire run, and the dashed line shows the spectrum for the wave forms truncated 
at t = 112. 5t D - 

FIG. 8. (a) The gravitational wave form rh + for Run 6 is shown for an observer located on 
the axis at = cf> = at distance r from the source, (b) The gravitational wave luminosity L/Lq, 
where Lq = c 5 /G. (c) The gravitational wave energy spectrum dE/df. The solid line shows the 
spectrum for the entire run, and the dashed line shows the spectrum for the wave forms truncated 
at t = 107.5t D - 

FIG. 9. Particle positions are shown projected onto the x — y plane for Run 7. Here, 
Mi = 1.4M , M 2 = 0.85Mi, R x = R 2 = 10 km, T = 2, and t D = 0.073 ms. The initial sep- 
aration ao = 4i?i. The stars are orbiting in the counter-clockwise direction. The vertical axis in 
each frame is y/R\ and the horizontal axis is x/R±. 

FIG. 10. (a) The gravitational wave form rh + for Run 7 is shown for an observer located on 
the axis at 6 = (f> = at distance r from the source, (b) The gravitational wave luminosity L/Lq, 
where Lq = c 5 /G. (c) The gravitational wave energy spectrum dE/df. The solid line shows the 
spectrum for the entire run, and the dashed line shows the spectrum for the wave forms truncated 
at t = 125t D - 
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FIG. 11. Particle positions are shown projected onto the x — y plane for Run 8. Here, 
Mi = 1.4M , M 2 = 0.85Mi, Ri = 10 km, R 2 = 9.7 km, T = 3, and t D = 0.073 ms. The 
initial separation ao = 4R±. The stars are orbiting in the counter-clockwise direction. The vertical 
axis in each frame is y/R\ and the horizontal axis is x/R±. 

FIG. 12. (a) The gravitational wave form rh + for Run 8 is shown for an observer located on 
the axis at 6 = <\> = at distance r from the source, (b) The gravitational wave luminosity L/Lq, 
where Lq = c 5 /G. (c) The gravitational wave energy spectrum dE/df. The solid line shows the 
spectrum for the entire run, and the dashed line shows the spectrum for the wave forms truncated 
at t = 112. 5t D - 

FIG. 13. Particle positions are shown projected onto the x — y plane for Run 9. Here, 
Mi = 1.4M , M 2 = 0.5Mi, R 1 = R 2 = 10 km, T = 2, and t D = 0.073 ms. The initial sepa- 
ration ao = 4Ri. The stars are orbiting in the counter-clockwise direction. The vertical axis in 
each frame is y/Ri and the horizontal axis is x/R\. 

FIG. 14. (a) The gravitational wave form rh + for Run 9 is shown for an observer located on 
the axis at 9 = <p = at distance r from the source, (b) The gravitational wave luminosity L/Lq, 
where Lq = c 5 /G. (c) The gravitational wave energy spectrum dE/df. The solid line shows the 
spectrum for the entire run, and the dashed line shows the spectrum for the wave forms truncated 
at t = 215t D - 

FIG. 15. Particle positions are shown projected onto the x — y plane for Run 10. Here, 
Mi = 1.4M©, M 2 = 0.5Mi, Ri = 10 km, R 2 = 8.7 km, T = 3, and t D = 0.073 ms. The ini- 
tial separation ao = 4.0i?i. The stars are orbiting in the counter-clockwise direction. The vertical 
axis in each frame is y/R\ and the horizontal axis is xjR\. At the end of the run, frame (i), the 
secondary has roughly 1/3 of its original mass. 
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FIG. 16. (a) The gravitational wave form rh + for Run 10 is shown for an observer located on 
the axis at = (f> = at distance r from the source, (b) The gravitational wave luminosity L/Lq, 
where Lq = c 5 /G. (c) The gravitational wave energy spectrum dE/df. The solid line shows the 
spectrum for the entire run, and the dashed line shows the spectrum for the wave forms truncated 
at t = 220t D - 
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